print('initializing general updating coefficients')
from  calculate_domain_size import nxp1,ny,nz,nx,nyp1,nzp1
from initialize_fdtd_parameters_and_arrays import *
from calculate_material_component_values import *
from define_problem_space_parameters import dx,dy,dz
import numpy as np
#     tau_k= 7e-10
psi_x = np.zeros((nx  , nyp1 , nzp1))
psi_y = np.zeros((nxp1, ny   , nzp1))
psi_z = np.zeros((nxp1, nyp1 , nz))
psi_x = eps_0 * deps_x
psi_y = eps_0 * deps_y
psi_z = eps_0 * deps_z
den_x = 2*tau_k_x+dt
den_y = 2*tau_k_y+dt
den_z = 2*tau_k_z+dt
ksi_x=2*dt*psi_x/den_x
ksi_y=2*dt*psi_y/den_y
ksi_z=2*dt*psi_z/den_z

Cexe = np.zeros((nx, nyp1, nzp1))
Cexhz = np.zeros((nx, nyp1, nzp1))
Cexhy = np.zeros((nx, nyp1, nzp1))

Ceye = np.zeros((nxp1, ny, nzp1))
Ceyhx = np.zeros((nxp1, ny, nzp1))
Ceyhz = np.zeros((nxp1, ny, nzp1))

Ceze = np.zeros((nxp1, nyp1, nz))
Cezhy = np.zeros((nxp1, nyp1, nz))
Cezhx = np.zeros((nxp1, nyp1, nz))

# General electric field updating coefficients
# Coeffiecients updating Ex
for ii in range(nx):
    for jj in range(nyp1):
        for kk in range(nzp1):
            if eps_r_x[ii, jj, kk] == 1:
                Cexe[ii, jj, kk] = (2 * eps_r_x[ii, jj, kk] * eps_0 - dt * sigma_e_x[ii, jj, kk]) / \
                                   (2 * eps_r_x[ii, jj, kk] * eps_0 + dt * sigma_e_x[ii, jj, kk])
                Cexhz[ii, jj, kk] = (2 * dt / dy) / \
                                    (2 * eps_r_x[ii, jj, kk] * eps_0 + dt * sigma_e_x[ii, jj, kk])
                Cexhy[ii, jj, kk] = -(2 * dt / dz) / \
                                    (2 * eps_r_x[ii, jj, kk] * eps_0 + dt * sigma_e_x[ii, jj, kk])
            else:
                Cexe[ii, jj, kk] = (2 * eps_r_x[ii, jj, kk] * eps_0 - dt * sigma_e_x[ii, jj, kk] + ksi_x[ii, jj, kk]) / \
                                   (2 * eps_r_x[ii, jj, kk] * eps_0 + dt * sigma_e_x[ii, jj, kk] + ksi_x[ii, jj, kk])
                Cexhz[ii, jj, kk] = (2 * dt / dy) / \
                                    (2 * eps_r_x[ii, jj, kk] * eps_0 + dt * sigma_e_x[ii, jj, kk] + ksi_x[ii, jj, kk])
                Cexhy[ii, jj, kk] = -(2 * dt / dz) / \
                                    (2 * eps_r_x[ii, jj, kk] * eps_0 + dt * sigma_e_x[ii, jj, kk] + ksi_x[ii, jj, kk])

# 计算 Ey 更新系数
for ii in range(nxp1):
    for jj in range(ny):
        for kk in range(nzp1):
            if eps_r_y[ii, jj, kk] != 1:
                Ceye[ii, jj, kk] = (2 * eps_r_y[ii, jj, kk] * eps_0 - dt * sigma_e_y[ii, jj, kk] + ksi_y[ii, jj, kk]) / \
                                   (2 * eps_r_y[ii, jj, kk] * eps_0 + dt * sigma_e_y[ii, jj, kk] + ksi_y[ii, jj, kk])
                Ceyhx[ii, jj, kk] = (2 * dt / dz) / \
                                    (2 * eps_r_y[ii, jj, kk] * eps_0 + dt * sigma_e_y[ii, jj, kk] + ksi_y[ii, jj, kk])
                Ceyhz[ii, jj, kk] = -(2 * dt / dx) / \
                                    (2 * eps_r_y[ii, jj, kk] * eps_0 + dt * sigma_e_y[ii, jj, kk] + ksi_y[ii, jj, kk])
            else:
                Ceye[ii, jj, kk] = (2 * eps_r_y[ii, jj, kk] * eps_0 - dt * sigma_e_y[ii, jj, kk]) / \
                                   (2 * eps_r_y[ii, jj, kk] * eps_0 + dt * sigma_e_y[ii, jj, kk])
                Ceyhx[ii, jj, kk] = (2 * dt / dz) / \
                                    (2 * eps_r_y[ii, jj, kk] * eps_0 + dt * sigma_e_y[ii, jj, kk])
                Ceyhz[ii, jj, kk] = -(2 * dt / dx) / \
                                    (2 * eps_r_y[ii, jj, kk] * eps_0 + dt * sigma_e_y[ii, jj, kk])

# 计算 Ez 更新系数
for ii in range(nxp1):
    for jj in range(nyp1):
        for kk in range(nz):
            if eps_r_z[ii, jj, kk] != 1:
                Ceze[ii, jj, kk] = (2 * eps_r_z[ii, jj, kk] * eps_0 - dt * sigma_e_z[ii, jj, kk] + ksi_z[ii, jj, kk]) / \
                                   (2 * eps_r_z[ii, jj, kk] * eps_0 + dt * sigma_e_z[ii, jj, kk] + ksi_z[ii, jj, kk])
                Cezhy[ii, jj, kk] = (2 * dt / dx) / \
                                    (2 * eps_r_z[ii, jj, kk] * eps_0 + dt * sigma_e_z[ii, jj, kk] + ksi_z[ii, jj, kk])
                Cezhx[ii, jj, kk] = -(2 * dt / dy) / \
                                    (2 * eps_r_z[ii, jj, kk] * eps_0 + dt * sigma_e_z[ii, jj, kk] + ksi_z[ii, jj, kk])
            else:
                Ceze[ii, jj, kk] = (2 * eps_r_z[ii, jj, kk] * eps_0 - dt * sigma_e_z[ii, jj, kk]) / \
                                   (2 * eps_r_z[ii, jj, kk] * eps_0 + dt * sigma_e_z[ii, jj, kk])
                Cezhy[ii, jj, kk] = (2 * dt / dx) / \
                                    (2 * eps_r_z[ii, jj, kk] * eps_0 + dt * sigma_e_z[ii, jj, kk])
                Cezhx[ii, jj, kk] = -(2 * dt / dy) / \
                                    (2 * eps_r_z[ii, jj, kk] * eps_0 + dt * sigma_e_z[ii, jj, kk])

# 计算磁场更新系数 (向量化操作)
Chxh = (2 * mu_r_x * mu_0 - dt * sigma_m_x) / (2 * mu_r_x * mu_0 + dt * sigma_m_x)
Chxez = -(2 * dt / dy) / (2 * mu_r_x * mu_0 + dt * sigma_m_x)
Chxey = (2 * dt / dz) / (2 * mu_r_x * mu_0 + dt * sigma_m_x)

Chyh = (2 * mu_r_y * mu_0 - dt * sigma_m_y) / (2 * mu_r_y * mu_0 + dt * sigma_m_y)
Chyex = -(2 * dt / dz) / (2 * mu_r_y * mu_0 + dt * sigma_m_y)
Chyez = (2 * dt / dx) / (2 * mu_r_y * mu_0 + dt * sigma_m_y)

Chzh = (2 * mu_r_z * mu_0 - dt * sigma_m_z) / (2 * mu_r_z * mu_0 + dt * sigma_m_z)
Chzey = -(2 * dt / dx) / (2 * mu_r_z * mu_0 + dt * sigma_m_z)
Chzex = (2 * dt / dy) / (2 * mu_r_z * mu_0 + dt * sigma_m_z)
def printInfo(fields,i,j,k):
    print("%.30f"%(fields[i-1,j-1,k-1]))
# Initialize coeffiecients for dispersive materials
from initialize_dispersive_coefficients import *

# Initialize coeffiecients for lumped element components
# initialize_voltage_source_updating_coefficients
# initialize_current_source_updating_coefficients
# initialize_resistor_updating_coefficients
# initialize_capacitor_updating_coefficients
# initialize_inductor_updating_coefficients
# initialize_diode_updating_coefficients
# initialize_thin_wire_updating_coefficients
# initialize_incident_field_updating_coefficients
